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We use a two-dimensional (2D) elastic free energy to calculate the effective interaction between 
two circular disks immersed in smectic-C films. For strong homeotropic anchoring, the distortion of 
the director field caused by the disks generates additional topological defects that induce an effective 
interaction between the disks. We use finite elements, with adaptive meshing, to minimize the 2D 
elastic free energy. The method is shown to be accurate and efficient for inhomogeneities on the 
length scales set by the disks and the defects, that differ by up to 3 orders of magnitude. We compute 
the effective interaction between two disk-defect pairs in a simple (linear) configuration. For large 
disk separations, D, the elastic free energy scales as ~ D~ 2 , confirming the dipolar character of 
the long-range effective interaction. For small D the energy exhibits a pronounced minimum. The 
lowest energy corresponds to a symmetrical configuration of the disk-deflect pairs, with the inner 
defect at the mid-point between the disks. The disks are separated by a distance that is twice the 
distance of the outer defect from the nearest disk. The latter is identical to the equilibrium distance 
of a defect nucleated by an isolated disk. 



I. INTRODUCTION 



Colloidal systems consist of a dispersive medium or 
solvent and a dispersed phase. The size of the dispersed 
domains or colloids is large compared to atoms but small 
compared to macroscopic lengths. Thus these systems 
exhibit new phenomena that are challenging theoretically 
while they are of practical importance in a variety of ap- 
plications, from paints and coatings to food and drugs 
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The colloidal state is metastable and its life-time is de- 
termined (largely) by the effective interactions between 
colloids. In many systems these interactions are tunable 
through the control of one or more external parameters. 
To prevent phase separation due to van der Waals or 
other attractive interactions, the colloids are coated or 
charged to induce steric or electrostatic repulsions |L|. 
Recently, a novel mechanism of colloidal interactions was 
reported: the elastic distortion of a liquid crystal host in 
a reverse nematic emulsion H, i.e., an emulsion of wa- 
ter droplets in a nematic liquid. This interaction arises 
from a competition between the aligning properties of 
the spherical droplet surfaces that favour e.g., a radial 
(homeotropic) nematic orientation and the bulk elastic- 
ity and boundary conditions that favour a uniform ne- 
matic. This effect is generic and it arises in a variety of 
anisotropic liquids for different isotropic inclusions. In 
addition, elastic mediated interactions were shown to be 
effective above the bulk nematic-isotropic temperature, 
for colloids wet by the nematic phase J3|. Finally, un- 
der appropriate conditions, inverted nematic emulsions 
exhibit colloidal self-organization, i.e., new phases where 
the colloids are ordered over macroscopic distances j|,[| . 

Understanding the phase behaviour of these complex 
systems begins by identifying the mechanisms and cal- 
culating the effective interactions between colloids. This 



is far from a simple task owing to the non-linear nature 
of the nematic elasticity and the existence of widely dif- 
ferent length scales: the bulk correlation length that sets 
the scale for variations of the microscopic degrees of free- 
dom, such as nematic order, biaxility etc., and the size of 
the colloids that sets the scale for variations of the elastic 
degrees of freedom, i.e., the nematic director field, n(r) 
[S|. Both microscopic (molecular dynamics simulations 
[7 1, tensor order-parameter free energy H) and macro- 
scopic (elastic free energy minimization using: trial func- 
tions H||, simulated anealing fll0|] , finite elements |ll]] 
and electrostatic analogies J12|) approaches have been 
used to calculate the orientational order of a nematic in 
the presence of a single colloid, in 2 and 3 dimensions. 

Theoretical |j{ll]] and experimental (||l3|,[l4| work on 
inverted nematic emulsions in 3D, reveals that for suffi- 
ciently strong homeotropic anchoring, spherical droplets 
induce topological defects, i.e., singularities in the ne- 
matic director field. There are several of these defects 
and their stability depends on parameters such as the size 
of the droplet, the anchoring strength and the boundary 
conditions. When the size of the droplet is large, a hy- 
perbolic point defect appears at a certain distance from 
the droplet, and this structure is known as the satellite 
configuration |Q-|ll| . For smaller droplets and/or weaker 
anchoring strengths, the minimum elastic free energy cor- 
responds to a ring disclination that surrounds the droplet 
on the equatorial plane, and this structure is known as 
the saturn-ring configuration ||-|l^]. The distortion of 
the director field due to the presence of more than one of 
these colloids induces multipolar long-range interactions 
between the colloids [g,|l5],[L6| . Depending on the sym- 
metry of the distortion, the dominant long-range inter- 
actions are dipolar [[15[ or quadrupolar 1 17 14| 
range, however, the colloids repel each other 
prevents their colapse. 
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A 2D system was also investigated by Pettey et al. |1J] 
who obtained analytic solutions for the smectic-C direc- 
tor around circular disks. For a disk at the center of 
a smectic-C layer, the solution that minimizes the elas- 
tic free energy exhibits a hyperbolic point defect, at a 
distance r° = \[2a from the center. The free energy of a 
system of multiple disk-defect pairs was also obtained, for 
pairs that are far apart. For defects aligned with the far 
field director, a dipole-dipole interaction results between 
the colloids Jf^] . These predictions were tested in a recent 
experiment by Cluzeau et al. |l8| . Thin films of smectic- 
C* were heated and the temperature was stabilized when 
cholesteric (N*) droplets started to nucleate. Cluzeau 
and collaborators observed that a hyperbolic point defect 
is nucleated with each cholesteric droplet, at the distance 
predicted by theory. In addition, they observed that the 
droplets interact through a long-range attraction and as- 
semble into chains, confirming the dipolar nature of the 
interaction. Finally, they measured a mean separation 
D = (1.3 ± 0.1)d a between nearest-neighbor droplets in 
chains, where d a is the droplet average diameter. 

Even in 2D exact solutions are not readily obtained. 
The analytic solutions of JT3] are restricted to one in- 
clusion or to multiple inclusions that are sufficiently far 
apart. Recently, Fukuda and Yokoyama solved the 2D 
problem using the tensor order parameter formalism and 
an efficient adaptive grid method to integrate numeri- 
cally the dynamic equation Q. One of the advantages 
of this method is that no special treatment of topologi- 
cal defects is required, by contrast with the director de- 
scription, where defects appear as singularities of the di- 
rector field H. For all the parameters used in m, the 
authors found that the equilibrium director configura- 
tion around a circular disk, exhibits quadrupolar sym- 
metry and corresponds to two defects with topological 
charge -1/2. This contrasts with the 3D result, where the 
quadrupolar defect was found to be stable over a limited 
range of parameters. It is also in apparent contradiction 
with the experiment on 2D smectic-C* films [lSy . The 
latter, however, do not exhibit mirror symmetry about 
reflection of the 2D director, excluding a distortion with 
quadrupolar symmetry |pi|,p|. 

In this paper we use finite elements and adaptive 
meshes to minimize the 2D elastic free energy functional. 
We focus on dipolar distortion fields as observed in the 
experiments frq ]. In particular, we study the interaction 
between two disks of equal radii, a, in a large smectic- 
C layer. The minimal free energy occurs when the two 
dipoles are aligned head to tail, in a direction parallel to 
the far field director. In this geometry we investigated a 
wide range of disk separations. We verified the validity 
of the electrostatic analogy at large distances and calcu- 
lated the short-range repulsive interaction between two 
disks. 



II. MODEL AND NUMERICAL METHOD 



The free energy of a nematic is invariant under uniform 
rotations of the whole sample and under the symmetry 
operations n — » — n and r — » r. In addition, n is a unit 
vector so that riiVjUi is zero. These considerations imply 
that the stiffness tensor has 2 independent bulk compo- 
nents in 2D. The distortions whose energy is measured 
by these constants are (1) splay, with non-zero V • n, and 
(2) bend, with non-zero V x n. The elastic free energy for 
the bulk nematic phase, is the sum of these terms p^,|S[: 



F = 



[K S {V ■ n) 2 + K B (V x n) 2 }dS 



(1) 



where K$ and Kb are the elastic constants associated 
with splay and bend distortions, and n = (cos 9, sin 9) is 
the 2D nematic director. The integral is over the region 
occupied by the 2D nematic film. 

Note that if K$ = Kb = K, the elastic free energy 
is proportional to ^- ■ V^n, V^n.,- (plus boundary terms) 
which is the form produced by simple models of the ne- 
matic state, such as those with pairwise intermolecular 
interactions depending only on the angle between the 
long molecular axes or the 2D XY-model [O]. 



F = y /[(V-nr :-(V> n) 2 ]dS 
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(2) 
(3) 



Line and boundary terms have been ignored in (y) which 
is justified in the strong anchoring regime (fixed orienta- 
tion at the boundaries) for systems with large line ten- 
sions (inclusions with fixed shape) since their contribu- 
tion to the free energy is constant. The one elastic con- 
stant approximation (0) is not strictly necessary, since 
the numerical method to be described below can be ap- 
plied to the Frank free energy (|l|); the choice of (j^) how- 
ever, allows a direct comparison with the analytical re- 
sults of Jl^] and provides a stringent test of our numerics. 

The geometry used in the numerical calculations is 
shown schematically in figure la. We consider two circu- 
lar disks of radius a separated by a distance D. Close to 
each disk we pin a defect at distances r<jj and rd 2 mea- 
sured from the center of each disk, respectively. The sys- 
tem has mirror symmetry with respect to x-axis. Thus, 
the center of the disks and the accompanying defects are 
located on the x-axis. The smaller circles around the 
defects are the defect cores (regions where the nematic 
order is destroyed [pf). They are used as cut-offs for the 
integral of the elastic free energy density (||). Lengths are 
measured in units of a and energy is measured in units 
of K. 
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FIG. 1. a : Schematic representation of the geometry with 
mirror symmetry about the a;-axis. b : Detail of the final 
adapted mesh, close to the defect core indicated by the solid 
line. The triangulation carefully respects the circular contour. 
c : The 2D elastic free energy as a function of the adaptive 
mesh iteration index. After 4 or 5 iterations, the relative error 
in the computed free energy is 10~ 4 . Figures b and c corre- 
spond to one isolated disk in a smectic C layer. Lengths are 
measured in units of the disk radius a and energy is measured 
in units of K. 

We fixed the nematic orientation, 9, at the physical 
boundaries. The far field director was taken parallel to 
the x axis (9 = 0) and homeotropic boundary conditions 
were imposed at the disk perimeters. The solution for the 
director field is found by minimizing the total free energy 
of the system, using standard numerical procedures p9[ . 
The major difficulty in the numerical problem stems from 
the large difference in the length scales set by the disk and 
the defect. It was overcome by using adaptive meshing 
techniques. With the help of a 2D mesh generator |2C[ ] 
(BL2D package), a first triangulation respecting the pre- 
defined physical boundaries, as well as the defect cores, is 
constructed. A final grid for the region close to a defect is 
shown in figure lb and was obtained after several adapta- 
tion loops. Note that in addition to the outer boundaries 
(figure la), the triangulation follows the perimeter of the 
defect core. Using standard numerical procedures the 
elastic free energy is minimized in a region that includes 
the defect core. The elastic free energy of the core re- 
gion is subtracted at the end of the calculation. The core 
energy was set to zero [||. 

Adapted meshes may be generated if the Hessian of 
the solution, d 2 9/dxdy, is known. Far from the disks, 9 
varies slowly and the triangles can be large. By contrast, 
close to the defects, the Hessian is large (it diverges at 
Td) and the triangles have to be made several orders of 



magnitude smaller. The Hessian, ?&, at vertex k is esti- 
mated [g3[ by using the following 'weak' definition ( J20| , 
p. 349) 



n_/k 



J OXi OXj 

Jv k dS ' 



(4) 



where the integral is over the region occupied by the 
smectic layer. v k is the piecewise linear hat function asso- 
ciated with vertex k (v k = 1 at vertex k and v k — else- 
where), x\ — x and X2 — y- To construct an anisotropic 
non-uniform mesh, we require a metric map that dis- 
tributes the interpolation error, i.e., the difference be- 
tween the interpolated function and the exact solution, 
in a uniform fashion fl20fl . The required map, A4 k ^, at ver- 
tex k is a two dimensional positive definite matrix given 
by 



M k = c a O 



|Ai| 
|A 2 | 



0" 



(5) 



where O is the orthogonal matrix that diagonalizes the 
Hessian H k , and Ai and A 2 the corresponding eigenval- 
ues. The constant cq controls the number of mesh trian- 
gles. To avoid problems caused by the divergence of the 
Hessian, the eigenvalues are bounded by 
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where ft„j„ and hf nax are the minimal and maximal mesh 
edge lengths. 

As a test of the numerics, we considered one single disk 
in a smectic-C film. We calculated the free energy as a 
function of the defect distance Td, for a dipole parallel 
to the far field director. The equilibrium free energy was 
found for a defect at r^ = 1.41 ±0.01, in good agreement 
with the analytical result for small defects (ra = v2) [I4J 
and in line with the experimental result for inclusions of 
cholesteric droplets in smectic-C* films (r<j = 1.4 ± 0.1) 

S 



III. RESULTS 

In the calculations, we used a rectangular integration 
region of size L x x L y = 50 x 20 (see figure la). With the 
adaptive mesh technique we have studied defects with 
cores as small as 10 -3 , with a minimal mesh edge length 
of 10~ 4 . During the adaptation loops this limit was 
reached close to and inside the cores, providing a good 
discretization of the core region (see figure lb). The final 
mesh configuration contained up to 10 4 vertices. After 4 
or 5 iterations of the adaptive mesh procedure, the rela- 
tive error of the free energy is of the order of 10 -4 (see 
figure lc). For a particular geometry, the central part of 
the final mesh is shown in figure 2a. In figure 2b, a gray 
scale was used to plot the corresponding equilibrium tilt 
angle, 9(x, y). White is 9 = and black 9 = tt. 
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FIG. 2. a : Central part of the final adapted mesh, for 
a disk separation corresponding to the minimum free energy 
(L x x L y = 50 x 20). b : In the same region, the equilibrium 
tilt angle, 8(x, y), is plotted using a gray scale. White is 6 — 
and black 9 = it. Lengths are in units of the particle radius 



Near the core (in a small region of size < 10~ 2 ) the 
elastic distortion field is determined by the defect alone. 
We have checked numerically that a change in the core 
size from £i to & leads to a free energy difference of 
7riTlog(^i/^2) IJJIl; independent of the position of the 
disks. 

For these systems (small defect cores) we found that 
the separation of the disks had no influence on the posi- 
tion of the outer defect, r^ : the minimum free energy 
was obtained for an outer defect at r°, the position of the 
defect for an isolated disk, i.e., independent of the pres- 
ence of the first disk. In what follows we set ra 2 = 1.41. 

The elastic free energy as a function of the position of 
the inner defect, r^ , is shown in figure 3a, for several disk 
separations. The lowest free energy was found for a disk 
separation D = 2.82, with an inner defect at r^ — 1.41. 
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FIG. 3. a : The 2D elastic free energy as a function of 
the position of the inner defect ra 1 , for several disk separa- 
tions. The lowest free energy corresponds to a disk separation 
D = 2.82 and an inner defect at rd 1 — 1.41. b : The position 
of the inner defect, r^ 1 , that minimizes the elastic free energy 
as a function of the disk separation, D. Lengths are in units 
of the disk radius a and energy in units of K. 

In figure 3b we plot the position of the inner defect, 
r^j , that minimizes the elastic free energy, as a function 
of the disk separation, D. We identify three regimes. For 
small disk separations (2 < D < 2.82) the position of the 
inner defect, r^, is symmetrical with respect to the disks 
and scales linearly with D, namely r^ = D/2, i.e., the 



defect is at the mid-point between the disks. In the in- 
termediate regime (2.82 < D < 6.5) the defect position 
varies non-monotonically with D, reaching a maximum, 
rd 1 = 1.52, at D = 3.5. Finally, for large disk separations 
(D > 6.5) the position of the inner defect is independent 
of D and is given by that of an isolated disk, r^ = 1.41, 
i.e., it is independent of the presence of the second disk. 

Note that the minimal free energy corresponds to a 
symmetrical configuration: the inner defect is at a dis- 
tance from the first disk identical to the distance of the 
outer defect from the second disk. 

Finally, we plot the elastic free energy as a function 
of the disk separation, D, in figure 4a. The free energy 
exhibits a pronounced minimum at D — 2.82±0.01, close 
to the value 2.6 ±0.2 observed experimentally by Cluzeau 
et al. JL8|, for chains of (polydispersed) circular droplets. 
Note that the experimental estimate of D, if taken with- 
out error bars, suggests that the inner defect is nucleated 
at a distance smaller than the distance corresponding to 
isolated droplet-defect pairs. This configuration, how- 
ever, is ruled out by our model (see figure 3a, open circles 
and squares). Within the (large) error bars we conclude 
that the experimental result is in line with the theoret- 
ical prediction: the lowest free energy corresponds to a 
symmetrical configuration, with the inner defect at a dis- 
tance from the disks that is identical to the distance of 
deffects in isolated disk-defect pairs. 

Admitedly, some effects were left out from our calcu- 
lation, the most important of which is arguably that of 
fluctuations. In fact, fluctuations were predicted to be 
important in 2D [12j but were reported to be ne glig ible 
in the rather thick films studied by Cluzeau et al. |18[. It 
is possible that many-body interactions may account for 
a reduction in disk separation for disks in chains, when 
compared to the equilibrium separation of a pair of disks. 
This could be checked within our model by calculating 
the free energy of three interacting disks. The experimen- 
tal results, however, suggest that this effect, if present, 
is of the order of the experimental error. A smaller disk 
separation, than that predicted by our model, may also 
result from surface contributions that were oversimplified 
in (0). In fact, deviations from strong homeotropic an- 
choring at the disks reduce the effective size of the inclu- 
sions, leading to a smaller equilibrium distance between 
disks [|§. 

The signature of the dipolar interaction for large disk 
separations may be seen in figure 4b. Here, the circles 
are the (calculated) force between disks, as a function of 
their separation. The force is clearly a power law that de- 
cays as D~ 3 . Note that the numerical results for the force 
are less accurate than those of the free energy, since they 
were obtained by numerical differentiation of the latter, 
figure 4a. 
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FIG. 4. a : The 2D elastic free energy as a function of the 
disk separation D. The free energy has a pronounced min- 
imum at D = 2.82 ± 0.01. b : Log-log plot of the effective 
attractive force between disks for large separations D (D > 5). 
Circles are the numerical values obtained by differentiation of 
the free energy shown in a. The line is D~ 3 . Lengths are in 
units of the disk radius a and energy is in units of K. 



IV. SUMMARY AND CONCLUSION 



We studied the interaction of two circular disks im- 
mersed in smectic-C films. For strong homeotropic an- 
choring the circular disks nucleate topological defects 
of the 2D director field. We implemented a finite ele- 
ments, with adaptive meshing, numerical method capa- 
ble of dealing with the variation of the director field over 
the widely different length scales set by the disks and the 
defects. The method allows the construction of well de- 
fined core regions around the defects and leads to very 
accurate results for the elastic free energy. 

We have verified the long-range dipolar character of 
the interaction between disk-defect pairs, predicted by 
Pettey et al. fl^j . In addition, we have calculated the 
interaction at small disk separations, and found that at 
static equilibrium, the distance between the centers of 
the disks is equal to D = 2.82 ± 0.01. 

Finally, we have determined the position of the inner 
defect, as a function of the separation between disks. For 
small disk separations the defect is located symmetrically 
at the mid-point between the disks. At intermediate sep- 
arations, this symmetry is broken and the position of the 
defect reaches a maximum r& x = 1.52 for a disk sepa- 
ration, D = 3.5. Finally, at large disk separations the 
defect relaxes slowly to the equilibrium position, r°, of 
an isolated disk. 

We intend to use the numerical method described in 
this paper to study the properties of other reverse ne- 
matic emulsions in 2 and 3D. In particular, we will con- 
sider the quadrupolar defects, recently studied by Fukuda 
and Yokoyama [g[ using a tensor order parameter formal- 
ism. Although these authors concluded that the dipolar 
configuration is dynamically unstable, and the quadrupo- 
lar configuration is the only one allowed in 2D nematics, 
the energetics of these two types of pinned defects has 
not been compared. 
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